clear all
set more off


global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"

cd $data

cap log close
log using $figures\L_file,replace t



***To get alpha
cap program drop get_alpha
program define get_alpha
	scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */
	gen xD=xbL_only*sqrt(1+s_omega^2+s_psi^2)+mu*female			
	ren a2 xA
	scalar rho_euA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)								/*We want the corr b/w epsilon and uA, not between uL and uA*/
	gen pDi=binormal(xA,xD,rho_euA)/normal(xA)
	gen m1i=normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,xD,rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,xD,rho_euA)^(-1)	
	gen m0i=-normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,-xD,-rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(1-normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,-xD,-rho_euA)^(-1)	
	scalar sigma=sqrt(1+s_csi^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)
	su part1 if female==1 & A==1
	scalar p1_1=r(mean)
	su part1 if female==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	su part2 if female==1 & A==1
	scalar p2_1=r(mean)
	su part2 if female==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	su part3 if female==1 & A==1
	scalar part3m=r(mean)
	gen rho_FX		=part1+part2
	gen lambda_FX	=part3*female
	gen F=female
	probit R female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9
	predict PI,pr
	gen psi=invnorm(PI)
	su psi 		if F==1 & R!=.
	scalar psi1=r(mean)
	su psi 		if F==0 & R!=.
	scalar psi0=r(mean)
	su rho 		if F==1 & R!=.
	scalar rho1=r(mean)
	su rho 		if F==0 & R!=.
	scalar rho0=r(mean)
	su lambda 	if F==1 & R!=.
	scalar lambda1=r(mean)
	scalar alpha1=((psi1-psi0)-(rho1-rho0))/lambda1		/*alpha D-in-D*/
	scalar list alpha1
end




u V_pars,clear
sort draw
merge draw using R_pars
drop if _merge!=3
drop _merge
sort draw
merge draw using A_pars
drop if _merge!=3
drop _merge

egen nonconv=rsum(cvg_VL cvg_ALR cvg_C cvg_0L cvg_0A cvg_0Aw_R cvg_partcorr cvg_chol)
keep if nonconv==8

set seed 110968
gen uu=uniform()
sort uu
keep in 1/200

drop uu draw

ren alpha1 alpha_old
ren alpha1_alt alpha1		/*This is the DD method*/
su alpha1
scalar se_alpha1=r(sd)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R alpha1,cov 
matrix V_alpha=r(C)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R,cov 
matrix V=r(C)
mat M3=vecdiag(V)'


***THIS IS EWMD (with the bootstrap var-cov matrix as weighting matrix)
u $data\moments_for_R,clear
mkmat M1
mata: mata clear
mata:
void eval0(todo, b, v, g, H)
 {
		  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
		  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
		  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
		  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
		  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
		  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
		  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
		  x8 =b[8]*b[4]+b[10]															
		  x9 =sqrt(b[8]^2+b[3]^2)														
		  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
		  x11=b[8]/sqrt(1+b[11]^2)
		  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
		  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
		  M1=st_matrix("M1")
		  M2=st_matrix("M3")
		  m1 =M1[1,1]
		  m2 =M1[2,1] 
		  m3 =M1[3,1]
		  m4 =M1[4,1]
		  m5 =M1[5,1] 
		  m6 =M1[6,1]
		  m7 =M1[7,1] 
		  m8 =M1[8,1]
		  m9 =M1[9,1]
		  m10=M1[10,1]
		  m11=M1[11,1]
		  m12=M1[12,1]
		  m13=M1[13,1]
	v=((x1-m1)^2/1)+((x2-m2)^2/1)+((x3-m3)^2/1)+((x4-m4)^2/1)+((x5-m5)^2/1)+((x6-m6)^2/1)+((x7-m7)^2/1)+((x8-m8)^2/1)+((x9-m9)^2/1)+((x10-m10)^2/1)+((x11-m11)^2/1)+((x12-m12)^2/1)+((x13-m13)^2/1)
 }
S = optimize_init()
optimize_init_evaluator(S, &eval0())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
optimize_init_which(S, "min")
optimize_init_which(S, "min")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
b = optimize(S)
b'
st_matrix("beta_ewmd",b')
end



***THIS IS DWMD (with the bootstrap diag of var-cov matrix as weighting matrix)
u $data\moments_for_R,clear
mkmat M1
mata: mata clear
mata:
void eval0(todo, b, v, g, H)
 {
		  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
		  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
		  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
		  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
		  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
		  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
		  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
		  x8 =b[8]*b[4]+b[10]															
		  x9 =sqrt(b[8]^2+b[3]^2)														
		  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
		  x11=b[8]/sqrt(1+b[11]^2)
		  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
		  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
		  M1=st_matrix("M1")
		  M2=st_matrix("M3")
		  m1 =M1[1,1]
		  m2 =M1[2,1] 
		  m3 =M1[3,1]
		  m4 =M1[4,1]
		  m5 =M1[5,1] 
		  m6 =M1[6,1]
		  m7 =M1[7,1] 
		  m8 =M1[8,1]
		  m9 =M1[9,1]
		  m10=M1[10,1]
		  m11=M1[11,1]
		  m12=M1[12,1]
		  m13=M1[13,1]
		  v1 =M2[1,1]
		  v2 =M2[2,1] 
		  v3 =M2[3,1]
		  v4 =M2[4,1]
		  v5 =M2[5,1]
		  v6 =M2[6,1]
		  v7 =M2[7,1] 
		  v8 =M2[8,1]
		  v9 =M2[9,1]
		  v10=M2[10,1]
		  v11=M2[11,1]
		  v12=M2[12,1]
		  v13=M2[13,1]
	v=((x1-m1)^2/v1)+((x2-m2)^2/v2)+((x3-m3)^2/v3)+((x4-m4)^2/v4)+((x5-m5)^2/v5)+((x6-m6)^2/v6)+((x7-m7)^2/v7)+((x8-m8)^2/v8)+((x9-m9)^2/v9)+((x10-m10)^2/v10)+((x11-m11)^2/v11)+((x12-m12)^2/v12)+((x13-m13)^2/v13)
 }
S = optimize_init()
optimize_init_evaluator(S, &eval0())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
optimize_init_which(S, "min")
optimize_init_which(S, "min")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
b = optimize(S)
b'
st_matrix("beta_dwmd",b')
end







scalar mu		=beta_ewmd[4,1]
scalar s_omega	=beta_ewmd[1,1]
scalar s_psi	=beta_ewmd[9,1]	
scalar s_v		=beta_ewmd[2,1]
scalar s_csi	=beta_ewmd[11,1]	
u data_to_est_alpha,clear
get_alpha
matrix beta=beta_ewmd\alpha1
svmat beta
keep beta
drop if beta==.
ren beta1 beta
gen str20 namepar=""
replace namepar="sigma_omega"      in 1
replace namepar="sigma_v"   	   in 2
replace namepar="sigma_phi"        in 3
replace namepar="pi"  		       in 4
replace namepar="gamma"  	       in 5
replace namepar="tau"   	       in 6
replace namepar="theta"  	       in 7
replace namepar="lamda"  	       in 8
replace namepar="sigma_psi"        in 9
replace namepar="psi"   	       in 10
replace namepar="sigma_csi"        in 11
replace namepar="alpha"        	   in 12
order namepar beta
save $data\struct_pars_estimates_ewmd,replace

scalar mu		=beta_dwmd[4,1]
scalar s_omega	=beta_dwmd[1,1]
scalar s_psi	=beta_dwmd[9,1]	
scalar s_v		=beta_dwmd[2,1]
scalar s_csi	=beta_dwmd[11,1]	
u data_to_est_alpha,clear
get_alpha
matrix beta=beta_dwmd\alpha1
svmat beta
keep beta
drop if beta==.
ren beta1 beta
gen str20 namepar=""
replace namepar="sigma_omega"      in 1
replace namepar="sigma_v"   	   in 2
replace namepar="sigma_phi"        in 3
replace namepar="pi"  		       in 4
replace namepar="gamma"  	       in 5
replace namepar="tau"   	       in 6
replace namepar="theta"  	       in 7
replace namepar="lamda"  	       in 8
replace namepar="sigma_psi"        in 9
replace namepar="psi"   	       in 10
replace namepar="sigma_csi"        in 11
replace namepar="alpha"        	   in 12
order namepar beta
save $data\struct_pars_estimates_dwmd,replace





u $data\struct_pars_estimates_ewmd,clear
su if namepar=="alpha"
scalar alpha=r(mean)
u $data\moments_for_R,clear
set obs 14
replace namevar="alpha" in 14
replace M1=alpha if namevar=="alpha"
replace M2=se_alpha1 if namevar=="alpha"
drop M2 namevar
svmat V_alpha
save $data\R_dataset_mom_se_ewmd,replace


u $data\struct_pars_estimates_dwmd,clear
su if namepar=="alpha"
scalar alpha=r(mean)
u $data\moments_for_R,clear
set obs 14
replace namevar="alpha" in 14
replace M1=alpha if namevar=="alpha"
replace M2=se_alpha1 if namevar=="alpha"
drop M2 namevar
svmat V_alpha
save $data\R_dataset_mom_se_dwmd,replace


/*
Now exit stata and run the R file $data\MD_EWMD.R and $data\MD_DWMD.R choose whether you want EWMD, DWMD -- This is to get s.e.'s and the GOF stat 
*/


cap log close
clear


